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ABSTRACT 

An expression is derived for calculating microcanonical-ensemble averages of the kinetic 
energies of atoms of different types in clusters isolated from the environment. This expres- 
sion is a natural generalization of the solution to the problem of hard spheres with different 
masses to a system with a many-particle interatomic interaction potential. The dynamics 
of a CgHg cubane is simulated numerically. The data on the numerical simulation confirm 
the validity of the results obtained. 



Statistical averages of physical quantities can be determined using ensembles of different 
types, such as microcanonical, canonical, large canonical, and other ensembles. In this 
case, the results obtained for macroscopic systems with numbers of particles ^ 1 differ 
by a negligibly small value of the order of 0(l/iV) [1]. Therefore, frequently (although 
not necessarily), the choice of a particular ensemble is determined by the convenience of 
calculations. In numerical simulation of the dynamics of finite systems, the number of 
particles and their total energy E (as well as the volume V occupied by these particles) 
often remain unchanged during the evolution, which corresponds to a microcanonical 
NVE ensemble [1,2]. In many models, the total momentum P is also an integral of 
motion (NVEP ensemble [3-5]. Prom the physical standpoint, this statement of the 
problem corresponds to systems that do not interact with the environment (for example, 
excited clusters in a gas phase [6-13]). 

If particles are classical and the energy of the interparticle interaction depends only 
on their relative coordinates (rather than on velocities), the kinetic energy per degree 
of freedom in a canonical ensemble (i.e., an ensemble of systems that are in thermal 
equilibrium with a reservoir at a temperature T) is equal to ^ksT, where ks is the 
Boltzmann constant. This is known as the equipartition theorem. Correspondingly, the 
kinetic energy of each particle (irrespective of its mass) is (£fej„)c = ^^bT, where i = 1 — N 
is the particle number (hereafter, it will be assumed that the system is three-dimensional) 
and the total kinetic energy is represented in the form 

{Ekin)c = ^kBTN, (1) 

where (...)c indicates the canonical-ensemble average. 

The relationship similar to expression (1) is valid for a microcanonical NVE ensemble 



[2]; that is, 
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{Ekin)NVE = -^ksT^nN, (2) 



where is the microcanonical temperature defined by the formula 
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Here 5* = /c^ In O is the entropy and VL[N, V, E) is the phase space volume proportional to 
the number of states with energies lower than or equal to E. Note that the definition of 
the entropy S through the logarithm of the density of states uj — dQ/dE leads to another 
relationship between the qmintities{E kin)NVE and T^, which transforms into expression 
(2) at A?" ^ oo (for more detail, see [2]). Despite the formal similarity of expressions 
(1) and (2), the equipartition in the NVE ensemble can be considered only "integrally" 
because it does not follow from relationship (2) that the kinetic energy (s^j^) nve averaged 
over the NVE ensemble is ^ksTm for each particle . 

If the total momentum of the system is P = 0, relationship (2) takes the form [3-5] 

{Ekin)NVEP = ^keT^iN - 1) (4) 

because the number of degrees of freedom decreases by three [note that the microcanonical 
temperatures in expressions (2) and (4) are different, since the additional condition 
P = results in a change in the phase space volume Jl] . When, in addition, the angular 
momentum is J = 0, we have 

{Ekin)NVEI>J = ^kBTrr^iN - 2), (5) 

because the number of degrees of freedom decreases by six. For ergodic systems, the 
microcanonical-ensemble average is equal to the average {■■■)t over the time of the evolution 
of one individual system. In this case, the microcanonical temperature coincides with 



the "dynamical temperature" that, at P = and J = 0, is determined from the formula 
[14] 

{Eun)t^\kBTa{N-2). (6) 

It should be noted that the use of periodic boundary conditions in numerical calculations 
(as a rule, in order to weaken the effect of finite sizes) is equivalent to the isolation of 
the system from the action of external forces. This leads to the conservation of the total 
momentum (but not the total angular momentum) [4]. 

Recently, Shirts et al. [15] succeeded in deriving the exact analytical relationships for 
the one-particle energy, momentum, and velocity distribution functions for a model system 
consisting of N hard spheres with E — const. In particular, it was demonstrated that 
the quantity {e\^^NVE — ^^bT^ in the NVE ensemble is identical for all atoms, whereas 
the equipartition theorem in the NVEY* ensemble with particles of different masses mj 
is violated; that is, 

(^fcm)wgp ^ Mtot — mg 

{^kin)NVEP Mtot — rUb 

where Mtot = Sili is the total mass of the system. In small systems, the ratio 
i^kin) NVEp/ i^kin) NVEP Can differ substantially from unity. This was confirmed in [15] 
using numerical calculations with periodic boundary conditions. Physically, the above 
effect is explained by the fact that the condition P = const differently affects the energy 
distribution of particles with different masses (for example, the maximum energy of one 
particle turns out to be dependent on the particle mass [15]). 

Shirts et al. [15] noted that their results can (possibly, with some corrections) appear 
to be valid for systems with realistic interatomic interaction potentials. In the present 
work, we first prove the equality of average kinetic energies of interacting particles in an 
NVE ensemble and derive formula (7) for interacting particles in an NVEP ensemble. 



Second, we demonstrate how this formula is changed when both the momentum and 
the angular momentum are conserved in the system. Third, we perform the numerical 
simulation of the dynamics of a CgHg cubane and show that the simulation results are in 
good agreement with the derived expression. Finally, it is demonstrated using the CgHg 
cubane as an example that the equipartition theorem holds true in the presence of the 
heat exchange with the environment. 

The average of the physical quantity ^(p, r) over an NVE ensemble can be represented 
in the form [2] 

[A) — — r„jv J — . f/c^ — 777 — v; — ' y°) 

lUlLi dpidriS{E-H{p,r)) 
where p = {pi} and r = {rj} are sets of momenta and coordinates of all particles, 
respectively, and if(p, r) = Z^iliPi/2mj + C/(r) is the Hamiltonian of the system. Let 
us determine the relation of the average {A)is[ve to the canonical-ensemble average A{T), 
which is given by the formula 

/ a-i dpidri exp (-^) A{p, r) 

A(T) = (A), = ' ^ t^y k^T .g. 

^ ^ ^ /nf.irfp.rfr,exp(-^) ^ ' 

By multiplying the integrands in the numerator and the denominator of relationship (9) 
by 1 = / dE6{E — H{p, r)), changing the order of integration over the energy E and the 
phase space, and using formula (8), we obtain 

^ JdEexp{-^)u;{E){A)NVE 
/rfEexp(-^)a;(E) ' 

where u;{E) is the denominator of formula (8), i.e., the density of states of the system. 

According to expression (10), the ratio between the averages A{T) and B{T) is written 

in the following form: 
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By setting A{p, r) = = p^/2ma and -B(p, r) — e^.^^ — pl/2mb and taking into account 
that, for the canonical ensemble, the equality {£kin)c = {^kin)c = ^^bT is true irrespective 
of the particle mass, we find from relation (11) that 

^kin)NVE 

The right-hand side of this equality does not depend on the temperature T. This is 
possible only when the equality {Skin)NVE — li^tin) nve is satisfied, where 7 is a constant 
independent of the energy E. It follows from formula (12) that 7 = 1, i. e., (£fej„)AryB — 
{s\ir^NVE at all energies E. Therefore, the kinetic energy average over the NVE ensemble 
is identical for all particles and, according to relationship (2), can be represented in the 

form {£kin)NVE = ^ksTm- 

For an NVEP ensemble, the average of the quantity A(p, r) is written in the form [3] 



^ / nil dp,dr,A{p, t)6{E - H{p, v))S (P - EIi Py 
IIll,dpdrAE-H{p,r))5{P-j:l,Pi) 

The canonical-ensemble average Ap(T) in the presence of the constraint on the total 

momentum is represented as follows: 

, _ _ J nL dp.dr. oxp {-SS2f) A(p, r)S (P - TSL, P.) 

" * " /n£.<iP.*.exp(-igfi)*(P-E£,P.) ■ ^ ' 

By using the same procedure as for the NVE ensemble, we obtain: 

^ JdEexp(-^).p(E)(^>...P 
/di5exp(-jfj)u-p(£) 

where u!p{E) is the denominator in formula (13). 

As a result, the ratio between the energies {£kin)c,F and {£kin)c,p is given by the formula 
(4.n)c,P _ /^^exp(-^)a;p(£;)(£^JwEP 



(4m)c,P JdEexp (-^) UJp{E){6li^)NVEP 
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In order to calculate the energy (£fej„)c,p, we change over in formula (14) to new Jacobi 
momentum variables [4, 16]: 

P^ = ^Ep^-#^P^+i> k<N-l, (17) 
-'^^fc+i i=i -'i^fc+i 

where = Sf=imi and Fn — J2iLiPi- Since the total kinetic energy is represented 
as Ekin = EiIiPf/2mi = E^=iPfc/2Aifc, where i^k = ruk+iMk/Mk+i ior k < N - 1 and 
/In — Mfot, and the Jacobian of the transformation from {pj} to {P^} is equal to unity 
[4, 16], from expression (14) at P = we derive the following relationship for the A^th 
particle: 



where we used the expression pjv = — Pat-i, which follows from formula (17). Since 
the Jacobi momenta depend explicitly on the atomic numbering and any particle can be 
chosen as the A^th particle, the following equality is satisfied: 

(4.n)c,P = IksT (l - ^) , i = 1 - TV . (19) 

Prom relationships (16) and (19), we have 

JdEe^ip (-^) u;p{E){eljNVEP ^ Mtpt - m„ 
JdEexp (-^) ujp{E){elJj,vEi> M,,t - m, " 

Since the right-hand side of this equality does not depend on the temperature T, we 

obtain formula (7) for the ratio (£fcj„)ArvEp/(^fcm) wbp- This means that the kinetic 

energy in the NVEP ensemble for an arbitrary interparticle interaction is distributed 

among particles in the same manner as in the model problem of hard spheres [15]. In this 

case, from relationship (4), we derive 

(4.n)wEP = hsT^ - ■ (21) 



The results obtained allow for the simple physical interpretation. Since the total mo- 
mentum in the NVE ensemble is not constant, the energy { ^^°^'^ )nve = ^ksT^, (where 
V is the velocity of the center of mass) corresponds to three degrees of freedom asso- 
ciated with the translational motion of the system as a whole. The contribution of the 
ith particle to the energy { '^*°2^'^ )nve is proportional to the particle mass and, hence, 
is ^ksTm^^- At P = (i. e., at V = 0) the average kinetic energy of the ith particle 
decreases by the above value and, therefore, is given by the expression |/cbTto(1 — ]^^)- 

By performing a similar analysis, it is possible to determine the average kinetic energy 
of the ith particle in the case of the additional constraint on the angular momentum J = 0. 
We assume that particles (atoms) form a cluster in which each particle executes vibrations 
with respect to its equilibrium position. At a small amplitude of these vibrations, to each 
of the three degrees of freedom associated with the rotation of the cluster as a whole in 
the NVE ensemble there approximately corresponds the energy { ^"^ ^'^ )nve — ^ksTm, 
where are principal moments of inertia of the cluster with respect to its center of 
mass in equilibrium and fin are components of the angular velocity {n — 1 — 3). The 
contribution of the ith particle to the energy ^" ) nve is proportional to its moment 
of inertia R and, hence, is ^ksT^^. As a result, at J = (i. e., at 17 = 0) the average 
kinetic energy of the ith particle decreases by Ylt,=i k^BTmj^ and (with allowance made 
for P = 0) is represented by the relationship 



Although relationship (22) is derived not rigorously but only for states close to equilibrium 
states (in which the moments of inertia can be approximately considered constant 
during the evolution), this relationship is in the excellent agreement with the molecular 
dynamics data. 




(22) 
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The CsHg cubane [17] was chosen for the numerical calculation of the distribution 
of the kinetic energy among atoms of different types. In this metastable cluster, the 
carbon atoms occupy cube vertices and the hydrogen atoms are located at extensions 
of the main diagonals (Fig. 1). We used the nonorthogonal tight-binding model for 
hydrocarbon systems [18], which was modified in our earlier work [19] in order to provide 
better agreement between theoretical and experimental values of the binding energies and 
interatomic distances in different molecules and clusters C^H^. For the bond lengths in 
the cubane, this model leads to Icc — 1-5696 A and Zch = 1.0823 A, which are close to 
the experimental bond lengths, to 1.571 A and 1.097 A, respectively. 

At the initial instant of time, random velocities and random displacements were assigned 
to each atom, so that the momentum and angular momentum of the whole system were 
equal to zero. Then, the classical Newton equations were numerically solved with the time 
step ^0 — 2.72 • 10~^^ c. The total cluster energy E — Epot + E^in remained constant in the 
course of simulation. The excitation energy E^^ equal to the difference between the total 
energy and the energy of the cubane in the equilibrium state was chosen not very high 
(in order to ensure the cluster lifetime long enough to collect sufficiently large statistics 
before decay of the cluster) and not very low (in order for the autocorrelation times 
determined by the anharmonicity of vibrations be as short as possible) . The determined 
optimum excitation energy E^x corresponds to a dynamical temperature ~ 1000 K 
[see relationship (6)]. In this case, (-E'j^.jjj)^ ~ 0.4925'g^ ^ 0.5Egx- This means that the 
anharmonicity effects, while weakly pronounced, occur in the system. 

Figure 2 shows two characteristic calculated dependences of the ratio between the evo- 
lution time-averaged kinetic energies (i.e., the kinetic energies averaged along the tra- 
jectory of cluster motion in the phase space) of the carbon and hydrogen subsystems 



on the number of molecular dynamics steps N. The statistical processing of the results 
leads to the ratios {Ekin,c)t/ {Ekin,B)t = 0.8153 ±0.0011 and 0.8161 ± 0.0008 at dynamical 
temperatures ~ 900 K and 1000 K, respectively, which, within the limits of error, coin- 
cide with the ratio {EI^j^^q) nvefj / {EUj^^b) nvepj = 0.8155, calculated from relationship 
(22) with due regard for the cubane parameters Mfot — 8mc + Smn, = 8/^ + 8/^^, 
^n/^n — + 75 tec''^' ~ 12mH. This good agreement between the computer 

simulation data and the result obtained from relationship (22) indicates that the system 
is ergodic (the average over the ensemble of systems is equal to the average over the 
trajectory of one system), on the one hand, and seems to be rather unexpected, on the 
other hand. Actually, relationship (22) was derived within the approximation that the 
moments of inertia of atoms are constant during the evolution (i.e., under the assumption 
that changes in the shape and sizes of the cluster are very small), whereas the relative 
changes in the bond lengths upon vibrations of the cubane can be rather large and reach 
~ 10% at Ta ~ 1000 K. It is interesting to note that, as follows from relationship (22), 
the constraint on the angular momentum can lead to a nonuniform distribution of the 
kinetic energy even in clusters composed of atoms with identical masses. 

We also simulated the dynamics of the CgHg cubane with allowance made for the 
cubane interaction with atoms of a buffer gas, which had a fixed temperature and played 
the role of a thermal reservoir. As was expected, the velocity distribution of hydrogen 
and carbon atoms after thermalization of the cluster became Maxwellian irrespective of 
the initial conditions. The difference between the kinetic energies of the hydrogen and 
carbon subsystems was absent. 

In conclusion, it should be noted that the nonuniformity of the kinetic energy distribu- 
tion in small thermally insulated clusters composed of atoms with different masses can be 

1 n 



significant and this factor should be taken into account when analyzing rapid processes 
associated with the decay of clusters and their interaction with each other. 
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Fig. 2a. Dependences of the ratio {Ekin,c)t/ {Ekin,}i)t between the evolution-time- 
averaged kinetic energies of the carbon and hydrogen subsystems in the CgHg cubane 
on the number of molecular dynamics steps N. The horizontal line corresponds to the 
value {Ekin,c) NVEPj / {Ekin,H) NVEPj — 0.8155 calculated from relationship (22). Dynami- 
cal temperature is = 887 K. 
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Fig. 2b. Dependences of the ratio {Ekin,c)t/ {Ekin,H)t between the evolution-time-averaged 
kinetic energies of the carbon and hydrogen subsystems in the CgHg cubane on the 
number of molecular dynamics steps N. The horizontal line corresponds to the value 
{Ekin,c)NVEPj/{Ekin,B)NVEPj — 0.8155 Calculated from relationship (22). Dynamical 
temperature is Td — 984 K. 



